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^ ! Abstract 

r~| ' The quantum instanton approximation is a type of quantum transition state theory that cal- 

' culates the chemical reaction rate using the reactive flux correlation function and its low order 

B '. 

D ■ derivatives at time zero. Here we present several path-integral estimators for the latter quantities, 

^ ; 

^ ! which characterize the initial decay profile of the flux correlation function. As with the internal 

c/3 ■ 

! energy or heat capacity calculation, different estimators yield different variances (and therefore 

'c/2 ! 

different convergence properties) in a Monte Carlo calculation. Here we obtain a virial(-type) esti- 

^ ■ 

mator by using a coordinate scaling procedure rather than integration by parts, which allows more 
y—i ' computational benefits. We also consider two different methods for treating the fiux operator, i.e., 

' local-path and global-path approaches, in which the latter achieves a smaller variance at the cost 

^ : 

• of using second-order potential derivatives. Numerical tests are performed for a one-dimensional 

i Eckart barrier and a model proton transfer reaction in a polar solvent, which illustrates the reduced 

' variance of the virial estimator over the corresponding thermodynamic estimator. 

o : 

>>: 
^ : 

Ph. 

> '■ 
X- 



1 



I. INTRODUCTION 



Developing an accurate and practical method for computing chemical reaction rates is 
one of the fundamental subjects of theoretical chemistry. In this regard the most successful 
approach is probably classical transition state theory (TST),i^ which has been applied 
widely to numerous reactions including biological systems such as enzyme catalysis.— The 
robustness of TST comes from its simplicity, i.e., the rate is determined solely from the 
free energy difference between the reactant and the activated complex. TST relies on the 
assumption of no "recrossing" trajectories through the dividing surface, which is usually 
valid at not too high temperature and for large dimensional systems. While successful in 
many cases, TST has the inherent deficiency that it accounts for no quantum effects, which 
needs to be addressed in order to treat low-temperature or light-atom transfer reactions. 
A conventional remedy to this problem is to add quantum corrections in a posteriori man- 
ner, e.g., by multiplying a tunneling factor that is computed along a prescribed tunneling 
path.^ Another strategy is to try to develop a quantum TST (QTST) by starting from the 
rigorous quantum rate expression and make some approximations for neglecting recrossing 
effects. Several such theories exist,— ^'''i^i^i^^i^'^i'^^i'^^i'^^i'^^ though there is in principle no unique 
formulation (in contrast to the classical case). 

The quantum instanton (Ql) approximatio n^^i^'^'i^^i^^i^^i^'^i^^i^^i^'^ is a recently developed 
theory for chemical reaction rates that is among the category of QTST. While the original 
derivation was based on the semiclassical "instanton" (periodic orbit in imaginary time) 
model,— the working rate expression can be understood roughly as the second-order cumu- 
lant (or Gaussian) approximation to the fiux-fiux correlation function. 



Cff(t) = Cff(O) + ^^(1(0)^2 + . . . ^ Cff(O) exp 



lC'ff(0),2 
2Cff(0) 



(1) 



for which the rate constant is given by 

^ 1- rdtc.it) . , ' (2) 

Lir Jo L^r \ ^ ^_Cg(0)/Cff(0) 

(see Sec. |n] for details). This approximation can be viewed as a quantum analog of the 
(classical) TST assumption in the sense that all possible oscillations in Cs{t) at later times 
(quantum re-crossing flux) are neglected. Test calculations show that this QI approximation 
gives a rate accurate to within ~10 % of the exact rate when the reaction is "direct", and 



also to within a factor of 2 even for cases in which significant recrossing is expected (e.g., 
the colhnear Cl+HCl reaction).— The computational merit of Eq. Q is that it is expressed 
wholly in terms of the Boltzmann operator, and thus it can be evaluated rigorously even 
for complex molecular systems using imaginary-time path integrals. A previous paper has 
presented such a scheme,— in which the factor Cf[{0)/Qr in Eq. Q is evaluated as the 
barrier height of a particular free energy surface, while the remaining factor is calculated as 
the statistical average of some estimating function over the transition-state path ensemble. 
This computational scheme has been applied successfully to several benchmark systems 
including gas-phase reactions such as CH4 + H — > CH3 + H2,^^ a model proton transfer 
reaction in a polar solvent,— and an isomerization reaction of pentadiene.— 

The purpose of this paper is to present an improved path-integral estimator for computing 
the QI rate. In particular, we focus on the factor Cfr(0)/Cfr(0) in Eq. Q which character- 
izes the initial decay profile of Cs{t). This quantity involves several different estimators 
because of the presence of the second time derivative. The estimator used in previous work 
was of "thermodynamic" type,-^ and its variance thus grows rapidly as a function of the 
number of path variables employed in the path integration. As with the internal energy or 
heat-capacity calculation , it should be possible to transform the thermo- 
dynamic estimator into a virial form in order to reduce the statistical error. In this paper we 
present such a scheme based on a coordinate scaling procedure, rather than integration by 
parts, which is based on the recent study by Predescu et alM^ and possesses the following 
computational benefits: (i) the transformation to a virial estimator is quite straightforward 
in contrast to integration by parts; (ii) one can use a finite-difference technique in order to 
avoid explicit calculation of potential derivatives in the virial estimator; and (iii) higher- 
order time derivatives of Cff(t) such as Cg'^^(O), Cf\{i), can also be generated with little 
modifications to the code, which can be used as input for more flexible approximations^i to 
the true Cs{t) than Eq. (0). 

The remainder of this paper is as follows: In Sec.|n]we summarize the working expression 
of the QI theory. In Sec. IIIII we first consider an "off-diagonal" average energy and derive 
its thermodynamic and virial estimators to describe the basic idea of coordinate scaling. In 
Sees. HVl and El we apply the scaling procedure to quantum time correlation functions in 
order to obtain a virial estimator for the reaction rate. In Sec. IVII we calculate the variance 
of the virial estimator for a one-dimensional Eckart barrier and a model proton transfer 
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reaction in a polar solvent. Sec. I VIII concludes. 



II. THE QUANTUM INSTANTON APPROXIMATION FOR CHEMICAL REAC- 
TION RATES 

The QI theory approximates the reaction rate as follows (see Ref. Q for the derivation 
motivated by semi classical considerations): 

Cff(0)v^ h 



where (0) is the zero time value of the fiux-fiux correlation function,^^ 

Cs{t) = tr L-PH/2p^-pH/2^iHt/Hp^-^Ht/H 

with F being the flux operator, 

{x — x^) + 5{x — x^)p], 
and AH represents a specific type of energy variance (i.e., AH^ = {H"^) — (H)^) 



(3) 



(4) 



F=-[H,hix-xi)] = — 
n 2m 



(5) 



AH' 







x^) 




He-P^l^ 








X 




{xt 


g-/3H/2 





n 2 



(6) 



In this paper we consider a one-dimensional system with the Hamiltonian H = /2m+V{x) 
for notational simplicity. In Eq. (jSJ, X"'- is the location of the dividing surface that separates 
the reactant and product regions. AH in Eq. (jH)) can be written more compactly as 

^'C'dd(O) 



where Cdd{t) is a "delta-delta" correlation function defined by 

Cdd(t) = tr \e'^^^^5{x - x^)e-^^^^e'^'^^5{x - x^)e-'^'/^^ 

Substituting Eq. ((7j) into Eq. Q gives 

Cff(O) 



(7) 



(8) 



k{T) 



Qr 



-C'dd(0)/Cdd(0) 



(9) 
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which has a formal resemblance to the Gaussian approximation to Cs{t) in Eq. (PJ. An 
extended version of the QI theory has also been proposed, which makes a log-augmented 
cumulant expansion of Cs(t) as follows,— 

2" 



Cff(t) ~Cff(0)exp<^6oln 



1 + 



hp) 



+ bit^ + ---bNt 



2N 



(10) 



where coefficients {bk} are determined by a matching procedure with the direct Taylor series 
expansion of Cfj(t). We note that the above approximation still falls among QTST because 
Cff(t) in Eq. (fTUj) is always positive and thus does not describe any recrossing effects (see 
Ref. |2l| for how this extension improves upon the QI rate). 



III. PATH-INTEGRAL ESTIMATORS FOR OFF-DIAGONAL AVERAGE EN- 
ERGY 

Before proceeding, it is useful first to present the coordinate scaling idea in its simplest 
form by considering an "off-diagonal" average energy defined by 

{Xb\e-P^\xa} op 

where ptaiP) = {xb\€-^^^\xa) , since this quantity serves as the basis for treating a time 
correlation function. Using the primitive approximation to the Boltzmann operator, 

^-(3H ^ L-eV/2^-ef^-eV/2Y (^2) 



with e = P/N, a. discretized path integral for is obtained as 

PbaiP) = j dXi--- j dXN-lWba{Xl, XN-I, P), (13) 



where 



WbaiXi, . . . , Xn-1; P) = ( o^^2^ J 



271} 



(14) 

with xo = Xa, xn = Xb, and Wk = 1/2 for k = 0,N and Wk = I otherwise. Differentiating 
Eq. ()13|) with respect to P gives a thermodynamic estimator for the energy, 

EUP) = {^T)ba (15) 



N ^ N 

^{xk - Xk-if + —^^WkVixk), (16) 



with 

2(3 2^2/^2 ""^-^^ ' 

where (■ ■ ■ )ba denotes an ensemble average over the weight function Wba{xi, . . . , x^-i, P)- 
This estimator has the well-known drawback that the statistical error grows with N due 
to cancellation of the first two terms in the right-hand side of Eq. As in the case of 

the internal energy or heat capacity, one can transform the above estimator into a virial 
form through integration by parts. Here instead we employ a coordinate scaling procedure 
that we find more useful.— >22iM To this end we first write the density matrix at a different 
temperature P', 

PbaiP') = dx[--- dx'^_,Wba{x',, . . . , /?'), (17) 



and then transform the integration variables {x'^} into a set of new variables {xk} according 
to 

x'k = xl + Jjixk- xl), (18) 
where x^ is the reference point given by 

k 

Xl = Xa + {Xb - Xa) — . (19) 

Using the following identity (or with the method described in Appendix), 

^ |:(4 - 4.)^ 4 - ^.-.f + - i) (20) 

one can rewrite Eq. (fT7|) as 

PbaiP') = j dXi--- j dXN~lWba{Xl, Xn-I, P)Rba{P'), (21) 

where 

RbaiP') = 7" ' exp -T7 [^(^'fc) - ^(^^)] • (22) 

{xb\e-P^ \Xa) [ ^ J 

Note that all the j3' dependence is now embedded in the Rba factor. Differentiating Eq. (|21|) 

with respect to /?' and taking the limit j3' ^ (3 gives a virial estimator 



EbaiP) = {ev)ba 



dp' 



(23) 

/3'=/5/ ba 
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with 



1 



m 



2j3 2h?l3' 



1 ^ 



k=0 



-{xk - xl)V'{xk) + V{xk) 



Alternatively, one may evaluate the virial estimator via finite difference asM 



26P 



(24) 



(25) 



in order to avoid explicit calculation of the potential derivatives. 



IV. VIRIAL ESTIMATOR FOR THE TIME DERIVATIVE OF CORRELATION 
FUNCTIONS 



With the scaling procedure above it is now straightforward to derive a virial estimator 
for the time derivative of correlation functions such as C*fr(0) and Cdd(O). We start with the 
following correlation function, 



C{t) = tr [e-^^/2^(a;)e-^^/2g^m/?i5(£)e-iH*/^" 



(26) 



where A and B are arbitrary position-dependent operators [note that C{t) becomes the 
delta-delta correlation function in Eq. (jH)) if we set A{x) = B{x) = 6{x — x^)]. For simplicity 
we work with the imaginary-time counterpart. 



C{X) = C{-ih\) = tr [e-(^/2+A)i/^(^)g-(/3/2-A)i?^^^) 



(27) 



with which the time derivative is given by {d/dt)^C{0) = {i/h)"'{d/dX)'^C{0). Discretizing 
the Boltzmann operators exp[—{(3/2 ± X)H] with P/2 time slices gives 



C(A) = J dxi - ■ ■ J dxpA{xo)B{xp/2)W{xi, . . . ,xp]X), 



where 



W{xi, ...,xp;X) 



mP 



2Tih\(3 + 2A) 



P/4 



mP 



2'Kh\(3-2X) 



P/A 



exp(— S') 



(28) 



(29) 



and 



S 



mP 



2h^{(3 + 2X) 



p/2 
k=l 



Xk Xk- 



P/2 

i)' + -^5^*fe(/? + 2A)\/(x,) 



k=0 



mP 
'2h\(5-2X) 



E 



[Xk Xk- 



k=P/2+l 



1 P 



(30) 



k=P/2 



with xo = xp and Wk = 1/2 for k = 0,P/2,P and Wk = I otherwise. Differentiating C{X) 
in Eq. ()28|) with respect to A and taking the A — hmit gives a thermodjTiamic estimator 
for C(0) (note that the first derivative vanishes by symmetry): 

(7(0) 



(31) 



where 



mP 



P/2 



P/2-1 



P-1 



k=l k=P/2+l 



P 



E- E Ki-') 

k=l k=P/2+l 



and 



2P 4mP^ , 



(3^ h?(3^ 



(32) 



(33) 



fc=i 



with ( ■ ■ ■ ) denoting an ensemble average over the weight function 
A{xo)B{xp/2)W{xi, . . . ,xp;0). This is the estimator that has been employed in pre- 
vious work .'^^i^?!^^!^? To transform it into virial form, we write C{X) in terms of temporary 
variables {x^}, 



C(A) = / dx[ ■ ■ ■ / dx'pA{xQ)B{x'p/2)W{x[, . . . , x'p] A) 



(34) 



and introduce a set of new variables {xk} as follows: 

'4 + ^^(2:^-4) iO<k<P/2) 

{P/2 <k<P) 
{k = 0,P/2,P) 

{0<k< P/2). 



^k~^ \/ ^ B i^k Xl) 



(35) 



with 



Xu — X I 



-k = XQ + ixp/2 - Xq) 



P/2 



The expression for C(A) then becomes 



C(A) = j dxi--- j dxpA{xQ)B{xp/2)W{xi,. . . ,xp;Q)R{X), 



(36) 



(37) 



where -R(A) = -Rkin-Rpot with 



-Rkin 



(xp|e-(^/2-^)^|a;p/2)(a;p/2|e-('5/2+^)^|a;o) 



{xp\e f^'^l'^\xp/2){xp/2\ 



-l3f/2 



Xo) 



(38) 



8 



and 

r . p/2 

Rpot = exp < --p JZ^I^fc + 2X)V{x',) - pV{xk)] (39) 



k=P/2 J 

Differentiating this expression for C{X) witli respect to A gives the desired virial estimator, 



with 




-1-IF} + G,) (41) 



^{xk - xl)V'{xk) + V{xk) 



(42) 



and 



16m 



p2 f^2p3 " ^P/^y ~]jpYl [^(^^ ~ xl)V'{xk) + (Xk - xlfV"{xk)] . (43) 

In practice we can avoid the calculation of first- and second-order potential derivatives by 
numerically differentiating R{X) as 

(7(0) 1 / R{6X) + R{-6X) -2R{0)\ 



c(o) ~ (sxy / ■ 



(44) 



V. TREATMENT OF THE FLUX OPERATOR 



A. Local-path approach 



Applying the above scheme to the flux-flux correlation function is somewhat tricky be- 
cause of the nonlocal character of the flux operator (i.e., a derivative operator). Different 
estimators arise depending on the route of the derivation, which in general exhibit different 
magnitudes of the variance. In previous wor k^''i'^^i^^i^^ we have employed a "local-path" es- 
timator, in which the flux operator was evaluated in terms of a few path variables near the 
dividing surface. This local estimator can be combined with the coordinate scaling proce- 
dure as follows. First we construct a discretized path integral for Cg{X) = Cf[{—ihX) as in 
Sec. IIVI in which the following matrix element appears: 



where {x'^} are temporary variables to be scaled later. Making the primitive approximation 
to e-(/3±2A)H/p evaluating the flux operator analytically via Eq. (0) gives 

ii-fi ~ y dx'.Six'o - a;*)t;o(A)(a;;|e-(^+2^)^/^|x[,)(x[,|e-('^-2^)^/^|x'„i), (46) 

where the velocity factor vo{X) is defined by 

with = 0. The effect of the flux operator is thus expressed in terms of only three path 
variables. Treating another flux operator in Cff(A) with the same method and performing 
the coordinate scaling precisely as in the preceding section gives 

C'ff(A) = J dxi--- J dxpW{xi, . . . ,xp;0)S{xo - x^)S{xp/2 - x^)R{\) (48) 

with 

R{X)=vo{X)vp/2{~X)R{X), (49) 

where W{xi, . . . ,xp;0) and -R(A) has the same definition as in Sec. IIVI Thus, the time 
derivative of Cs{t) can be obtained as 

(7ff(0) ^ _]_/ RjSX) +R{-6X) -2R{0) 
Cdd(O) " ~¥ \ (5A)^ 

Similarly, virial estimators for higher time derivatives, d"'Cg{0)/dt"'{n = 4,6,...), can be 
generated using an appropriate finite-difference formula of higher order.— 




B. Global-path approach 



One can also devise an alternate "global-path" estimator by first performing the coordi- 
nate scaling and then applying the flux operator (i.e., in an opposite order to the preceding 
section). To be specific, we insert the coordinate representation of the flux operator. 



2mi 



dx 



\x){x\ + \x){x\ 



(51) 
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with |x) = d\x)/dx into the imaginary-time flux correlation function as 



(/3/2- 




\xb){xb\e~ 


^(/3/2+A)H| 


Xa) 


{/3/2- 


-X)H 


\xb){xb\e~ 


^(/3/2+A)i?| 




-(/3/2- 


-\)H 


\xb){xb\e~ 


-(/3/2+A)H 


\Xa) 



where 



and 



(52) 



which can be written more compactly as 

Cs{X) = j dXa j dxb6{xa- x^)6{xb- x^) 

where an operator representing the "square" of the flux operator is given by 

.„„ Urn i^^^ ?1 (54) 

\2mJ x^^xax^^xt [dx+dxj ox~dXf^ dx+dx^^ dx^dxj ) 

Next we use a generalized scaling relation of the form (see Appendix): 

^x'b\e-'''^K) = /^^i---/ dxN.,Wba{xi,...,xr,-i;(3)R*{(3'), (55) 



; / I —/3'f\ r \ f 1 



I -I P' 
Xk = x,^ + J —{xk — xl), (57) 

4 = < + (58) 



Other quantities such as Wba and are deflned the same as in Sec. IIIII We note that the 
end-points (x^, xj,) are included in the coordinate transformation in addition to j3'. Applying 
the above relation to the density matrix elements in Eq. ()53p with N = P/2 and appropriate 
choice of end-points gives 



C'ff(A) 



j dxi--- j dxpS{xo - x^)5{xp/2 - x^)W{xi, ...,xp; 0)J^^R*{x^, x^, A), (59) 
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where R*{x^, x^,\)= Rt^Rfot with 



-"-kin 





e 


-(/3/2-A)T 




e 


-i(3/2+X)T 




(a;p|e-/3'^/2|. 


T^P/2){xp/2 


e-PT/2 





(60) 



and 



P/2 



-'^pot 



k=0 



{P-2\)V{x^')-pV{xk) 



k=P/2 



X 



X, 



Xi 



l(3±2X 



xZ + [xl - X, 



[Xk 

k 



^k)y 



Xh + ix„ - X 



k-P/2 
P/2 ' 



(61) 



(62a) 
(62b) 
(62c) 



The kth coordinate in Eq. (j62j) with plus and minus signs are defined for < k < P/2 and 
P/2 < k < P, respectively. The time derivative of Cfj(t) can be obtained by differentiating 
the factor jF^_R*(x^, x^, A) with respect to A, where the JF^ operator is applied analytically 
using up to second-order potential derivatives.^^ The latter operation is costly but often 
not too demanding because JF^ involves only the coordinates that define the (generalized) 
reaction coordinate, e.g., only a few Cartesian coordinates that describe the reacting atoms. 



VI. NUMERICAL TESTS 



We now apply the above estimators to a one- dimensional system with the Eckart potential 
barrier, 

V{x) = Vosech'^{x/a), (63) 

where Vq = 0.425 eV, a = 0.734 au, and the mass is 1060 au, which corresponds roughly to 
the H+H2 reaction. Table 1 lists the statistical error of Cf^/C^d and C^^/C^d = 2,4,6) 
obtained with 1 million path samples (note that the time arguments are always t = and 
are omitted hereafter). Three estimators are compared: the thermodynamic estimator, the 
local-path virial estimator in Sec. IV At and the global-path virial estimator in Sec. IV Bl The 
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latter two differ only in the treatment of the flux operator. The number of path variables 
used was P = 8 for T =1000 K and P =40 for T =200 K, which have a discretization error 
of ~2 % of the exact (P — oo) value. The dividing surface was always set at the top of the 
barrier with X"'- = in Eq. (0). The reader is referred to Refs. andl34 on how these time 
derivatives can be used to improve the approximate rates. 

We see from Table 1 that the virial estimators always exhibit a smaller statistical error 
than the thermodynamic estimator, as expected. Between the two virial estimators, the 
global-path version has a smaller variance than the local one by using more information 
on the entire path. The exceedingly small errors of the global-path estimator (< 0.1 %) 
at T = 1000 K are somewhat fortuitous, because at this temperature the system is close 
to the free-particle limit and the global-path estimator becomes exact for a free particle 
irrespective of the number of path variables.— This situation does not occur for the local- 
path virial estimator, where the velocity factor in Eq. (jTTj) must be averaged even for a 
free particle to give the correct result. Another important fact is that the variance of the 
virial estimators is nearly independent of the order of time derivatives in contrast to the 
thermodynamic estimator, which agrees qualitatively with the previous study by Predescu 
for the same system using a Fourier-like path integral.— 

Figure 1 plots the statistical error of Ci^V^dd and ds^/Cdd at T = 200 K function of 
the number of path variables P. The variance of the thermodynamic estimator grows rapidly 
with P, and the growth rate is especially large for Cq \ The local-path virial estimator also 
exhibits an increasing variance, which is caused by the appearance of P in the numerator 
of the velocity factor in Eq. (jTfj) . The global-path virial estimator, on the other hand, 
has a nearly constant variance regardless of the value of P, thus facilitating the systematic 
convergence to the P ^ oo limit. 

Next we apply the present method to a model proton transfer reaction^i in a polar solvent, 
AH -\- B +HB^ , where A, H, and B represent a hydrogen-bonding complex dissolved 

in liquid methyl chloride at T = 250 K. The details of the model is given in Ref. 0. Here 
we quantize only the proton degree of freedom with P = 40 and use the path integral 
Monte Carlo (MC) scheme described in Ref. [2^ Figure 2 shows the convergence of C^d/Cdd, 
Cjp^ /Cdd = 2,4) as a function of MC cycles. In all cases the virial estimators outperform 
the thermodynamic estimators in convergence rate. In particular, the convergence of Cdd 
is very rapid when using the virial estimator, which is beneficial in calculating the QI rate 
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in Eq. Q- On the other hand, the statistical error becomes larger for /C^d, and it 
was difficult to converge with 2 million path samples for n > 6. This is in contrast to the 
one-dimensional Eckart barrier studied above, where the variance of the virial estimator 
was nearly independent of the order of time derivatives. Apart from differences in the 
dimensionality of the system, the variance may be increased by stiff potential walls in the 
solute potential (defined with Morse-like functions),-^ because the virial estimator for 
depends implicitly on the higher-order potential derivatives. For example, the local-path and 
global-path virial estimators for C^^^ depend on 7th- and 8th-order potential derivatives, 
although the numerical calculation by finite difference needs only the 1st- and 2nd-order 
derivatives of the potential. It is not clear at present to what extent this behavior is common 
for other potentials (including polynomial potentials). Nevertheless, the fast convergence of 
2nd time derivatives even for the present stiff potential is very encouraging when considering 
future applications of the QI theory to more complex chemical reactions in condensed phases. 

VII. CONCLUDING REMARKS 

Our main purpose in this paper has been to show how a virial estimator for the time 
derivative of correlation functions can be obtained straightforwardly via a coordinate scaling 
procedure, and that the resulting estimator has an expected smaller variance than the ther- 
modynamic estimator. We have also presented two methods for treating the flux operator, 
i.e., local-path and global-path approaches, in which the latter has a smaller variance. The 
second time derivative of C^dit) and Cfj(t) are clearly the most important quantities for the 
QI rate in Eq. Q or in Eq. An open problem is how to best utilize the higher-order 
derivatives in order to improve the accuracy of approximate rates. While some progress has 
been made in this direction,— more studies would be useful if we consider the availability 
of C^\o) at least for systems with well-behaved potentials. 
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APPENDIX: USING THE FEYNMAN-KAC FORMULA 

The scaled expression (j21|) can also be obtained as follows. Utilizing the integration 
variables {i/k} defined by 

x'k = + J ^^Vk, (A.l) 
V m 

one can transform Eq. ()17j) as follows, 

PbaiP') = (x,|e-^'^|x,)Eexp|-^^^i;fcVK)| (A.2) 

I k=o J 

with 

Jdyi---J dyN-i exp { -f Ef=i(l/fe - Vk-if }(■■■) 

E(---) = —r , (A.3) 

Jdyi---J dyN-i exp |-f Ef=i(2/fc - Vk-iY^ 

which becomes the Feynman-Kac formula in the N oo limit with {y^} representing the 
standard Brownian bridge. Rewriting the above equation as 

pUP') = (a:f,|e-^^|x,)Eexp|-|:^w;,V^(a;,)|i?,,(/?0, (A.4) 



fc=0 



where Rba{(3') is defined by Eq. (j221) and 



Xk = + \ yk, (A.5) 

V m 

and changing integration variables from {yk} to {xk} results in Eq. ()21|) . Combining 
Eqs. flA.l|) and ()A.5|) gives the coordinate transformation in Eq. (fTHj) . Similar procedures 
can be used to obtain a generahzed expression in Eq. (fS^ . 
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FIG. 1: Relative statistical error (%) of the thermodynamic, local-path virial, and global-path 
virial estimators for (a) C^^ / Cdd and (b) C^^ / C^d computed for the one-dimensional Eckart barrier 
at T = 200 K. P is the number of path variables. 

FIG. 2: Statistical convergence of the thermodynamic, local-path virial, and global-path virial 
estimators for (a) C^J/C^d, (b) Cg^ /Cdd, and (c) /Cdd computed for a model proton transfer 
reaction in a polar solvent. 
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Fig. 1 (a) 
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Fig. 1 (b) 
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Fig. 2 (a) 
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Fig. 2 (b) 
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Fig. 2 (c) 
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TABLE I: Relative statistical error (%) of the thermodynamic, local-path virial, and global-path 
virial estimators for the one-dimensional Eckart barrier. One million paths are sampled with 8 and 
40 path variables for T = 1000 and 200 K, respectively. 

thermodynamic virial (local-path) virial (global-path) 
T = 1000 K 

Cdd/Cdd 0.5 0.024 0.024 

ci^V<^dd 1-2 0.34 0.012 

CffVCdd 2.3 0.33 0.014 

cP/Cm 2.6 0.33 0.015 

T = 200 K 

c£VCdd 1.1 0.27 0.27 

cI^/Caa 3.3 1.4 0.41 

C^'^VCdd 8.1 1.7 0.44 

C^^VCdd 25 2.0 0.51 
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